
--[[===========  太阳天文数据计算  ==========
]]


-- 引入必要的库文件(lua编写), 内部库不需要require
require 'sys'

local GMST =0
local Solar_A,Solar_H= 0,0

local Solar ={
	JD		= 0,
	A		= 0,
	H		= 0,
	Adeg	= 0,
	Hdeg	= 0,
	Dip_Meridian = 0    -- 太阳对子午线的倾角

		-- 函数 --
	-- InclineRad( A, H)	--根据高度和方位角计算倾角（单位：弧度）
	-- InclineDeg( A, H)
	-- Date2JulianDay(year, month, day, hour, minutes, seconds) -- 计算UT时间对应的儒略日,返回： 低精度的JD儒略日， JD整数，JD小数
	-- Position( JD1, JD2, longitude, latitude) --根据儒略日计算太阳的方位角和高度角
}

--[[--------------------------------------------------------
--		根据高度和方位角计算倾角（单位：弧度）
------------------------------------------------------------
	function incline( A,  H) 
	A： 方位角度，北0度，东90度,南180°
	H： 方位角度
	返回： 倾角角度,左0度，垂直90度，右180度
----------------------------------------------------------]]--
function Solar.InclineRad( A, H)
	return  math.atan2( math.sin(A)/math.cos(H), math.sin(A) )
end

-- 方位角和高度角计算倾角（单位： 角度）
function Solar.InclineDeg( A, H)
	local hh = math.rad(H)		-- 将角度转换位弧度
	local aa = math.rad(A)
	aa = math.atan2( math.sin(hh)/math.cos(hh), math.sin(aa) )
	aa = math.deg(aa)
	while aa < -90 do aa = aa +360 end
	while aa > 270 do aa = aa -360 end

	return aa
end


-------------------------------------------------------------
-- 计算UT时间对应的儒略日,返回： 单精度的JD儒略日， JD整数，JD小数
-------------------------------------------------------------
function Solar.Date2JulianDay(year, month, day, hour, minutes, seconds)
    local Y = year
    local M = month
    --1，2月，相当于去年的第13月和14月
    if month <=2   then
        Y = Y-1
        M = M+12
	end
	
    -- //将日,时,分换算到天(D) = day + hour/24.0 + minutes/(24.0*60.0) + seconds/(24.0*60*60)
    local D2 = hour/24.0 + minutes/(24.0*60.0) + seconds/(24.0*60*60)
	--local D = day + D2
    local B = 2 - math.floor(Y/100)+ math.floor(Y/400)

    -- JD = math.floor(365.25*(Y+4716)) + math.floor(30.600001*(M+1))+ day + B + D2 - 1524.5
	--计算JD小数部分
	local JD,JD1,JD2=0,0,0
	if D2 >=0.5 then
		JD2 = D2-0.5
		JD1 = math.floor(365.25*(Y+4716)) + math.floor(30.6*(M+1))+ day + B -1524
	else 
		JD2 = 1+D2-0.5	--向整数借1天
		JD1 = math.floor(365.25*(Y+4716)) + math.floor(30.600001*(M+1))+ day + B -1524 -1
	end
	Solar.JD = JD1+JD2		--单精度
	Solar.JD1 = JD1
	Solar.JD2 = JD2

    return Solar.JD, Solar.JD1, Solar.JD2
end

-----------------------------azimuth angle 方位角---------------------------
--function SolarPos( JulianDay, longitude, latitude)  根据儒略日计算太阳的方位角和高度角
-- 返回： 方位角和高度角（角度），方位角和高度角（弧度）
--------------------------------------------------------
function Solar.Position( JD1, JD2, longitude, latitude)
	if not JD1 or not JD2 or not longitude or not latitude then 
		log.info("Solar.Position( JD1, JD2, longitude, latitude)","参数异常！")
		return 0,0,0,0 
	end

	local d=0

	-- T 是J2000.0 起算的儒略世纪数,T=(JD-2451545)/36525
	local T = (JD1 - 2451545)/36525 + JD2/36525 ;
	--log.info(" T=", T)

	--太阳几何平黄经（λ ） Lo = 280°.46645 + 36000°.76983*T + 0°.0003032*T^2 （Date平分点起算）
	d = 280.46645 + 36000.76983*T + 0.0003032*T*T
	d = math.fmod(d,360)
	local L0 = math.rad(d)

	--太阳平近点角M = 357°.52910 + 35999°.05030*T- 0°.0001559*T^2-0°.00000048*T^3 
	d = 357.5291 + 35999.0503*T - 0.0001559*T*T-0.00000048*T*T*T;
	d = math.fmod(d,360);
	local M = math.rad(d);


  	--地球轨道离心率e = 0.016708617 - 0.000042037*T;	// - 0.0000001236*T^2 
	local e = 0.016708617 - 0.000042037*T - 0.0000001236*T*T;

	--太阳中间方程 C = (1°.914600 - 0°.004817*T -0°.000014*T*T) * sin(M)+(0°.019993 - 0°.000101*T) * sin(2M)+ 0°.000290*sin(3M) 
	d = math.sin(M);
	d = (1.9146 - 0.004817*T -0.000014*T*T) * math.sin(M)+(0.019993 - 0.000101*T) * math.sin(2*M)+ 0.00029*math.sin(3*M);
	d = math.fmod(d,360);
	local C = math.rad(d);

	--太阳的真黄经ZHJ,弧度: Θ = L0 + C
	local ZHJ = L0 + C;

	--真近点角  v = M + C
	local v = M + C;

	--日地距离 R =1.000001018 (1-e^2) / (1+e*cos(v)), 后面忽略了,没有用到R
	local R = 1.000001018*(1-e*e)/ (1+e*math.cos(v));

	--章动修正及光行差修正OMG:  Ω = 125°.04 - 1934°.136*T
	d = 125.04 - 1934.136*T;
	-- d = fmod(d,360);
	local OMG = math.rad(d);

	--    视黄经SHJ: λ = Θ - 0°.00569 -0°.00478*sin(Ω), 后面忽略了,没有用到R
	local SHJ = ZHJ- math.rad(0.00569) - math.rad(0.00478)*math.sin(OMG);

	-- 平黄赤交角E0=εo = 23°26'21".448 - 46".8150*T - 0".00059*T^2 +0".001813*T^3 (21.2)式中T 是J2000.0 起算的
	d = 23.4392911-0.0013004167*T-0.000000163889*T*T+0.000000503611*T*T*T;
	local E0 =  math.rad(d);

	--月亮的平黄经Ly =L′ = 218.3164591 + 481267.88134236T - 0.0013268T*T+ 0.0000019T*T*T
	d = 218.3164591 + 481267.88134236*T - 0.0013268*T*T+ 0.0000019*T*T*T;
	d = math.fmod(d,360);
	local Ly =  math.rad(d);

	-- 交角章动 De,弧度   
	d = (9.2* math.cos(OMG)+0.57* math.cos(OMG*L0)+0.1* math.cos(2*Ly)-0.09* math.cos(2*OMG));
	d = d/3600;
	local De = math.rad(d);

	--真黄赤交角EE =Eo+dE(De)
	local EE = E0 + De

	--地心赤纬DXCW (δ)  :  sinδ=sinε sinΘ
	d = math.sin(EE) * math.sin(ZHJ)
	local DXCW = math.asin(d)

	--地心赤经 DXCJ, (α) :  tanα = cosε sinΘ / cosΘ
	d = math.cos(EE) * math.sin(ZHJ)/math.cos(ZHJ)

	if d<0 then
    d = math.atan2( math.cos(EE) * math.sin(ZHJ), math.cos(ZHJ)) + math.pi*2
	else
    d = math.atan2(math.cos(EE) * math.sin(ZHJ), math.cos(ZHJ))
	end
	local DXCJ = d;
	--log.debug("DXCJ=",DXCJ)

	--[[格林尼治平恒星时 GMST, θo = 280.46061837 + 360.98564736629*(JD-2451545.0)+ 0.000387933*T^2 - T^3/38710000
		计算过程: 360.98564736629*(JD-2451545.0) + 280.46061837
			= 0.98564736629*(JD-2451545.0) + 360*(JD-2451545.0	) + 280.46061837  --忽略360整数
			= 0.98564736629*(JD-2451545.0) + 360*JD2 + 280.46061837
			=  0.98564736629*(JD1+JD2-2451545.0) + 360*JD2 + 280.46061837
			=  0.98564736629*(JD1-2451545) + 0.98564736629*JD2 + 360*JD2 + 280.46061837
			]]-- 
	-- d = 280.46061837 + 360.98564736629*(JulianDay -2451545)+ 0.000387933*T*T - T*T*T/38710000;
	local d1 =  math.fmod( 0.98564736629*(JD1-2451545), 360)
	local d2 =  math.fmod( 0.98564736629*JD2 + 360*JD2 + 280.46061837, 360 )
	GMST = math.fmod(d1+d2, 360)
	--log.debug("GMST=",GMST)

	--时角(弧度) H =θo - L - α  = 格林尼治平恒星时GMST - (L 观者站经度) -DXCJ地心赤经
	--H = 本地时角，从南向西测量。如果θ 是本地恒星时，θo 是格林尼治恒星时，L 是观者站经度（从格林尼治向西为正，东为负），那么本地时角计算如下：H = θ - α 或 H =θo - L – α(α是赤经)
	d = math.deg(DXCJ) ;
	d = math.fmod(d+3600,360);
	d =  GMST - longitude - d;
	local H =  math.rad(d);


	--//太阳方位角A : tan(A)=sin(H)/(cos(H)*sin(φ)-tan(δ)*cos(φ)) 其中: φ= 观测者(站)纬度,北半球为正,南半球为负 ;  视赤纬SCW ( δ)
	d = math.sin(H);
	d = math.sin(math.rad(latitude));
	d = math.sin(DXCW) / math.cos(DXCW);

	d = (math.cos(H) * math.sin(math.rad(latitude)) - math.sin(DXCW)/math.cos(DXCW)*math.cos(math.rad(latitude)))
	d = math.atan2( math.sin(H), d); --这个方位数据，南为0，东+90度，西-90度
  	--#                      //d = PI/2 - d;          这个方位数据，东为0，南+90度，西180度
	d = math.pi + d;          -- 这个方位数据，北为0度，东为90，南+180度，西270度
	Solar_A = d

	-- 太阳高度角h :  sin(h)=sin(φ)*sin(δ)+cos(φ)*cos(δ)*cos(H)
	d = math.sin(math.rad(latitude))*math.sin(DXCW)+math.cos(math.rad(latitude))*math.cos(DXCW)*math.cos(H);
	d = math.asin(d);
	Solar_H = d
	--返回方位角和高度角（弧度） 
	-- log.debug("天文","计算结果 儒略日:", JD1,"+", JD2, " 经度=",longitude," 纬度=",latitude," 方位角=",  math.deg(Solar_A),"高度角=",math.deg(Solar_H))
	Solar.A = Solar_A
	Solar.H = Solar_H
	Solar.A_deg =  math.deg(Solar_A)
	Solar.H_deg =  math.deg(Solar_H)

	return Solar.A_deg, Solar.H_deg, Solar.A, Solar.H
	--return math.deg(Solar_A), math.deg(Solar_H),Solar_A , Solar_H
end

return Solar

